R Notebook

This is my script to try and analyse data from scratch - DATASET 2 -> MERSCOPE

library(Seurat)

Warning message:
replacing previous import ‘RcppML::nmf’ by ‘NMF::nmf’ when loading ‘CARD’ 

library(ggplot2)
library(scCustomize)
library(readr)
library(pheatmap)
library(matrixStats)
library(spdep)
library(geojsonR)

Read in DAta

# Add the image slot to the object - we've left out the segmentation data as it took ages to load 
coords <- CreateFOV(coords = list(centroids = CreateCentroids(data$centroids[Cells(seurat),])),
                    type = c("centroids"),
                    assay = "MERSCOPE")
# Add the image slot to the object - we've left out the segmentation data as it took ages to load 
coords <- CreateFOV(coords = list(centroids = CreateCentroids(data$centroids[Cells(seurat),])),
                    type = c("centroids"),
                    assay = "MERSCOPE")
seurat[["COLON"]] <- coords  # adding image has a warning check the image to see what it looks like 
ImageDimPlot(seurat,fov='COLON')
ImageDimPlot(seurat,fov='COLON')

# Mark the blank probes that are detected in each cell - don't think we have the other codeword categories that are present 
seurat[["Negative.Control.Codeword"]] <- CreateAssayObject(counts =data$transcripts[grepl('Blank',rownames(data$transcripts)),])
saveRDS(seurat,file='/project/shared/spatial_data_camp/HACKATHON/the_bigger_fish/cg_seurat_obj')

ImageFeaturePlot(seurat, "nCount_MERSCOPE") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

explore the

Look at number of genes per cell

ImageFeaturePlot(seurat, "nFeature_MERSCOPE", axes = TRUE) + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

The image is big and takes ages so lets make a subet to get out code correct on and work through

Crop function only crops the image not subset the whole thing

cropped <- Crop(seurat[[“COLON”]], x = c(2000, 5000), y = c(4000, 7000), coords = “plot”) seurat[[“ROI1”]] <- cropped

seurat_subset <- seurat[,seurat$global_X < 5000 & seurat$global_X >2000 & seurat$global_Y < 7000 & seurat$global_Y > 4000]

Warning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objectsWarning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objects

seurat_subset <- seurat[,seurat$global_X < 5000 & seurat$global_X >2000 & seurat$global_Y < 7000 & seurat$global_Y > 4000]

Warning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objectsWarning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objects

ImageDimPlot(seurat_subset, axes=T)

ImageFeaturePlot(seurat_subset, "nCount_MERSCOPE") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

ImageFeaturePlot(seurat_subset, "nFeature_MERSCOPE") + scale_fill_viridis_c() 

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

Interesting this is quite different to the Xenium data - there are a large number of cells with few features - maybe we can plot these somehow

cOuld filter and also chec

quantile(seurat_subset$nFeature_MERSCOPE, c(0.01, 0.1, 0.5, 0.9, 0.99))

 1% 10% 50% 90% 99% 
  3  53 130 197 254 

# Can use MERSCOPE Volume incase 
ImageFeaturePlot(seurat_subset, "volume") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

Compared to visium there looks like cells start with a larger volume - some prefiltering of small cells

ggplot(seurat_subset[[]], aes(nCount_MERSCOPE, volume)) + geom_point() 

Do minimum filter based on what can cluster on -

seurat_subset[["SIZE_FILTER_LARGE"]] <- seurat_subset$volume < quantile(seurat_subset$volume, .99)
ImageDimPlot(seurat_subset, group.by="SIZE_FILTER_LARGE")

seurat_subset[["SIZE_FILTER_SMALL"]] <- seurat_subset$volume > quantile(seurat_subset$volume, .01)
ImageDimPlot(seurat_subset, group.by="SIZE_FILTER_SMALL")

p1 <- VlnPlot(seurat_subset, "nFeature_MERSCOPE", group.by = "SIZE_FILTER_SMALL", pt.size = .1, alpha = .5) + labs(title="Small Cell Filter")

Warning: Default search for "data" layer in "MERSCOPE" assay yielded no results; utilizing "counts" layer instead.

p2 <- VlnPlot(seurat_subset, "nFeature_MERSCOPE", group.by = "SIZE_FILTER_LARGE", pt.size = .1, alpha = .5)+ labs(title="Large Cell Filter")

Warning: Default search for "data" layer in "MERSCOPE" assay yielded no results; utilizing "counts" layer instead.

p1 + p2

seurat_subset[["SIZE_FILTER_SMALL"]] <- seurat_subset$volume > quantile(seurat_subset$volume, .1)
ImageDimPlot(seurat_subset, group.by="SIZE_FILTER_SMALL")

seurat$TRANSCRIPT_FILTER <- seurat$nCount_MERSCOPE >= 15
ImageDimPlot(seurat, group.by="TRANSCRIPT_FILTER")

Some areas have lots of low transcripts becuase haven’t cleared or lots of mucus so large areas discarded

seurat_subset$TRANSCRIPT_FILTER <- seurat_subset$nCount_MERSCOPE >= 15
ImageDimPlot(seurat_subset, group.by="TRANSCRIPT_FILTER")

ImageFeaturePlot(seurat_subset, "nCount_Negative.Control.Codeword") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

Check the code words and see if its specific ones that are - and check the actual assingment of the negative control genes - can look at control probe versus real gene ratio in each cell and if high proportion from control probes only bin the gene

For now lets just carry on

Might want to adjust this and probe filter - see above either set as level or as blank / real transcript ratio

seurat_subset$PROBE_FILTER <- seurat_subset$nCount_Negative.Control.Codeword == 0
ImageDimPlot(seurat_subset, group.by="PROBE_FILTER")

Explore QC Metrics

seurat_filtered <- subset(seurat_subset, PROBE_FILTER & SIZE_FILTER_LARGE & SIZE_FILTER_SMALL & TRANSCRIPT_FILTER)

Warning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objectsWarning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objects

ImageDimPlot(seurat_filtered)

seurat_filtered <- SCTransform(seurat_filtered, assay = "MERSCOPE", clip.range = c(-10, 10))

Running SCTransform on assay: MERSCOPE
Running SCTransform on layer: counts
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Variance stabilizing transformation of count matrix of size 500 by 14374
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 493 genes, 5000 cells
Found 33 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 500 genes
Computing corrected count matrix for 500 genes
Calculating gene attributes
Wall clock passed: Time difference of 4.960257 secs
Determine variable features
Centering data matrix

  |                                                                                                              
  |                                                                                                        |   0%
  |                                                                                                              
  |========================================================================================================| 100%
Getting residuals for block 1(of 3) for counts dataset
Getting residuals for block 2(of 3) for counts dataset
Getting residuals for block 3(of 3) for counts dataset
Centering data matrix

  |                                                                                                              
  |                                                                                                        |   0%
  |                                                                                                              
  |========================================================================================================| 100%
Finished calculating residuals for counts
Set default assay to SCT

seurat_filtered <- RunPCA(seurat_filtered)

PC_ 1 
Positive:  EPCAM, EPHB3, SOX9, CDH1, RGMB, AXIN2, PKM, CTNNB1, VEGFA, MYC 
       LGR5, IDH1, CDCA7, ASCL2, PROX1, MUC1, ERBB2, EPHB4, JUN, ERBB3 
       PCNA, CDK4, MKI67, HDAC1, XBP1, NOTCH1, MET, MCM2, LGALS9, CA9 
Negative:  COL1A1, ACTA2, COL4A1, PDGFRB, COL5A1, MMP11, FN1, SOCS3, CD248, ETS1 
       DUSP1, IL1B, ITGA1, ENG, PLVAP, PECAM1, ANGPT2, VCAM1, NFKBIA, COL11A1 
       PDGFRA, CAV1, CSF3R, NRP1, ELN, TGFB1, ADAMTS4, S100A9, MYH11, NDUFA4L2 
PC_ 2 
Positive:  COL1A1, ACTA2, FN1, PDGFRB, COL4A1, MMP11, COL5A1, CD248, EPHB3, CTNNB1 
       RGMB, ITGA1, ITGB1, SOX9, EPCAM, SMOC2, ETS1, LGR5, CDH1, ANGPT2 
       AXIN2, COL11A1, ELN, VCAM1, NOTCH1, PDGFRA, IDH1, CDCA7, TGFBI, CDK4 
Negative:  IL1B, SOCS3, FOS, NFKBIA, CXCL8, CSF3R, S100A9, PTGS2, ITGAX, CXCL5 
       SOD2, JUNB, DUSP1, SPP1, PROK2, FFAR2, CCL3, ICAM1, FCGR2A, CD14 
       FCGR3A, CEBPB, SERPINA1, VEGFA, TLR2, LYZ, ITGB2, HLA-DRA, IL6R, CD83 
PC_ 3 
Positive:  PECAM1, PLVAP, VWF, ENG, KDR, MMRN2, PDGFB, ETS1, FLT4, CDH5 
       INSR, CLEC14A, COL4A1, CXCR4, ITGA5, NRP1, SERPINE1, FLT1, PREX2, WWTR1 
       ITGB2, HLA-DRA, ANGPT2, TNFRSF4, CD40, ADAMTS4, IL3RA, CLDN5, PTPRC, CD3E 
Negative:  COL1A1, MMP11, COL5A1, PDGFRB, IL1B, ACTA2, PDGFRA, COL11A1, WNT5A, CXCL8 
       PTGS2, FN1, ELN, LRP1, COL6A3, FOS, FBLN1, EGR1, PDPN, TNC 
       SOCS3, VCAM1, IRS1, JUNB, SOD2, CSF1, NFKBIA, CXCL5, SMO, FZD7 
PC_ 4 
Positive:  COL4A1, PLVAP, IL1B, FOS, VWF, FN1, SOCS3, KDR, MMRN2, PECAM1 
       ENG, ANGPT2, PDGFB, INSR, CDH5, PTGS2, ADAMTS4, CSF3R, NFKBIA, FLT4 
       PDGFRB, CLEC14A, SERPINE1, CD248, CXCL8, S100A9, FLT1, PGF, FFAR2, NDUFA4L2 
Negative:  HLA-DRA, ITGB2, PTPRC, TRAC, CXCR4, CYBB, CSF1R, CD3E, CD14, CD2 
       C1QC, HLA-DRB1, HLA-DQA1, CD4, CIITA, ZAP70, SPP1, CD3D, IL2RB, HLA-DMA 
       HLA-DPB1, LYZ, ITK, MRC1, CD28, HAVCR2, ITGAM, GATA3, FOXP3, CCR4 
PC_ 5 
Positive:  VEGFA, EPHA2, CEACAM1, LAMC2, ATF3, LAMB3, PPARD, KIT, EGR1, FOS 
       CA9, LRP1, MUC1, TCF7L2, CXCL1, DDIT3, SLC26A3, CDKN1A, JUNB, MET 
       MUC2, SPP1, PLOD2, NFKB2, STAT6, LDHA, BCL2L1, AMOTL2, BMP1, PKM 
Negative:  FN1, SMOC2, EPHB3, LGR5, RGMB, SOCS3, CTNNB1, CDCA7, PCNA, CSF3R 
       CDK4, S100A9, PROX1, IDH1, SOX9, IL1B, MCM6, BIRC5, SOD2, MCM2 
       CCNB1, ASCL2, TP53, SOX2, PLK1, MKI67, HLA-DRA, NOTCH1, FOXM1, HIF1A 

ElbowPlot(seurat_filtered, 50)

NA
NA

Filter poor quality cells

Cluster the data

PC_Plotting(seurat_filtered, dim_number = 1)

FeaturePlot(seurat_filtered, "COL1A1", reduction = "pca") + scale_color_viridis_c()

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

ImageFeaturePlot(seurat_filtered, "PC_1") + scale_fill_viridis_c()

Warning: No FOV associated with assay 'SCT', using global default FOVScale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

ImageFeaturePlot(seurat_filtered, "COL1A1", size=.5) + scale_fill_viridis_c()

Warning: No FOV associated with assay 'SCT', using global default FOVScale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

seurat_filtered <- RunUMAP(seurat_filtered, dims = 1:20)

16:01:44 UMAP embedding parameters a = 0.9922 b = 1.112
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by ‘BiocGenerics’
16:01:44 Read 14374 rows and found 20 numeric columns
16:01:44 Using Annoy for neighbor search, n_neighbors = 30
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by ‘BiocGenerics’
16:01:44 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
16:01:47 Writing NN index file to temp file /project/.tmpRsessions/RtmpQlyPKi/file19fb8d30796043
16:01:47 Searching Annoy index using 1 thread, search_k = 3000
16:01:54 Annoy recall = 100%
16:01:56 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
16:01:59 Initializing from normalized Laplacian + noise (using RSpectra)
16:02:00 Commencing optimization for 200 epochs, with 627536 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
16:02:10 Optimization finished

seurat_filtered <- FindNeighbors(seurat_filtered, reduction = "pca", dims = 1:20)

Computing nearest neighbor graph
Computing SNN

seurat_filtered <- FindClusters(seurat_filtered, resolution = 0.7)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 14374
Number of edges: 496411

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8493
Number of communities: 13
Elapsed time: 3 seconds

DimPlot(seurat_filtered, label=T, repel=T)

NA
NA

ImageDimPlot(seurat_filtered, size=.6)

Warning: No FOV associated with assay 'SCT', using global default FOV

Transfer the labels onto Xenium from sc RNA reference dataset - what works well and what does not

markers <- FindMarkers(seurat_filtered, ident.1="0", max.cells.per.ident=500)
head(markers)

FeaturePlot(seurat_filtered, "CD3E", label=T, repel=T)+ scale_color_viridis_c(direction=-1)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

FeaturePlot(seurat_filtered, "MS4A1", label=T, repel=T)+  scale_color_viridis_c(direction=-1)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

FeaturePlot(seurat_filtered, "KIT", label=T, repel=T)+ scale_color_viridis_c(direction=-1)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

markers <- FindAllMarkers(seurat_filtered, max.cells.per.ident = 500)

Calculating cluster 0
Calculating cluster 1
Calculating cluster 2
Calculating cluster 3
Calculating cluster 4
Calculating cluster 5
Calculating cluster 6
Calculating cluster 7
Calculating cluster 8
Calculating cluster 9
Calculating cluster 10
Calculating cluster 11
Calculating cluster 12

head(markers)

top <- Extract_Top_Markers(markers, num_genes = 5, named_vector = FALSE, make_unique = TRUE)
top

 [1] "SLC26A3"  "VEGFA"    "CEACAM1"  "KIT"      "ATF3"     "GZMB"     "PCNA"     "CDK4"     "BAX"      "NRAS"    
[11] "IL1B"     "CSF3R"    "PTGS2"    "SOCS3"    "CXCL8"    "SMOC2"    "LGR5"     "EPHA4"    "RET"      "EPHB3"   
[21] "PDGFRB"   "CD248"    "COL5A1"   "COL1A1"   "ACTA2"    "PLK1"     "CCNB1"    "MKI67"    "BIRC5"    "AURKB"   
[31] "BRCA1"    "MCM2"     "E2F1"     "MCM6"     "CCNE1"    "FN1"      "SOX2"     "PTGDR2"   "CLCA1"    "MUC2"    
[41] "SERPINA1" "CTLA4"    "SPP1"     "MMP12"    "MMP9"     "CXCR4"    "IL4I1"    "CSF1R"    "HLA-DRA"  "CIITA"   
[51] "HLA-DQA1" "HLA-DRB1" "PLVAP"    "VWF"      "PECAM1"   "MMRN2"    "FLT4"     "CD3E"     "TRAC"     "CD2"     
[61] "FOXP3"    "TRBC1"   

Clustered_DotPlot(seurat_filtered, features = top, k=18)

[[1]]

[[2]]

ImageDimPlot(seurat_filtered, cells=WhichCells(seurat_filtered, expression = seurat_clusters %in% c(0, 5)))

Warning: No FOV associated with assay 'SCT', using global default FOV

ImageDimPlot(seurat_filtered, axes = T)

Warning: No FOV associated with assay 'SCT', using global default FOV

ImageDimPlot(seurat_filtered, fov="ROI1", border.color = "black" )

how much of it is limitations of the panel

ref <- readRDS("/project/shared/spatial_data_camp/datasets/SINGLE_CELL_REFERENCES/COLON_HC_5K_CELLS.RDS")
DimPlot(ref)

ref <- SCTransform(ref, residual.features =rownames(seurat))

Running SCTransform on assay: RNA
Computing residuals for the 550 specified features
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 19008 by 5725
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 5000 cells
Found 80 outliers - those will be ignored in fitting/regularization step

Skip calculation of full residual matrix
Calculating gene attributes
Wall clock passed: Time difference of 14.36988 secs
Determine variable features
Setting min_variance based on median UMI:  0.16
Calculating residuals of type pearson for 477 genes

  |                                                                                                              
  |                                                                                                        |   0%
  |                                                                                                              
  |====================================================                                                    |  50%
  |                                                                                                              
  |========================================================================================================| 100%

Centering data matrix

  |                                                                                                              
  |                                                                                                        |   0%
  |                                                                                                              
  |========================================================================================================| 100%
Place corrected count matrix in counts slot
Set default assay to SCT

ref <- RunPCA(ref)

PC_ 1 
Positive:  SOCS3, CCL2, FOS, ICAM1, JUNB, JUN, EGR1, IFITM1, CXCL12, CCL8 
       COL1A1, FBLN1, CDKN1A, CAV1, PLVAP, MMP2, BST2, SOD2, ENG, DUSP1 
       FN1, MYC, HLA-DRB1, GPX3, CLEC14A, PDGFRA, CSF1, PECAM1, CXCL1, ACTA2 
Negative:  SLC26A3, EPCAM, MZB1, CEACAM1, PLA2G2A, DERL3, PKIB, CCL5, MUC4, CD79A 
       CDH1, CD3E, LGALS9, FABP2, CD177, TNFRSF17, CD27, LAMC2, LAMB3, SOX9 
       MUC1, PTPRC, CCL3, EPHA2, CDCA7, CA7, CCL28, IDH1, CCL4, ERBB3 
PC_ 2 
Positive:  HLA-DRB1, HLA-DRA, HLA-DPA1, HLA-DPB1, PLVAP, HLA-DQA1, PECAM1, HLA-DMA, CD83, CAV1 
       FLT1, CLDN5, CLEC14A, ENG, THBD, TGFBR2, ICAM2, CXCR4, VWF, IL3RA 
       KLF2, CD79A, HLA-B, MMRN2, CDH5, PTPRC, IFITM1, CCL4, KDR, SNAI1 
Negative:  CCL2, FOS, CCL8, EGR1, CXCL1, FBLN1, COL1A1, JUN, FN1, MMP2 
       ATF3, MYC, JUNB, SOCS3, CCL11, GPX3, PDGFRA, LRP1, VCAM1, CXCL2 
       EPCAM, CEBPB, MAFB, ACTA2, CSF1, COL5A1, PTGS2, COL6A3, PLA2G2A, SLC26A3 
PC_ 3 
Positive:  PLVAP, CAV1, PECAM1, FLT1, ENG, CLEC14A, CLDN5, SLC26A3, THBD, TGFBR2 
       VWF, IL3RA, EPCAM, ICAM2, IFITM1, KLF2, MMRN2, CDH5, CEACAM1, KDR 
       ANGPT2, SNAI1, PKIB, PLA2G2A, NRP1, WWTR1, CX3CL1, CXCL12, COL4A1, ADAMTS4 
Negative:  HLA-DQA1, CD83, CXCR4, HLA-DPB1, CD79A, CCL4, CCL3, CCL2, HLA-DRA, HLA-DPA1 
       MZB1, CCL5, PTPRC, HLA-DRB1, CCL8, DERL3, LYZ, C1QC, CD3E, SOD2 
       HLA-DMA, CCR7, MS4A1, CD68, IL1B, PTGS2, FBLN1, CD27, CXCL8, ITGB2 
PC_ 4 
Positive:  MZB1, DERL3, CCL5, TNFRSF17, CD27, CD3E, CCL3, CCL4, CD79A, IRF4 
       CAV1, POU2AF1, CCR10, HLA-B, TNFRSF18, IFITM1, KLRB1, ICAM2, PLVAP, TGFB1 
       TNFRSF4, IL23A, PTPRC, CDKN1A, TNFRSF13B, PECAM1, CD3D, CTSW, NKG7, KLF2 
Negative:  CXCL2, HLA-DRA, HLA-DQA1, HLA-DPB1, HLA-DRB1, CD83, HLA-DPA1, SLC26A3, CXCL1, EPCAM 
       CXCL8, LYZ, HLA-DMA, PLA2G2A, CEACAM1, PKIB, CD68, ATF3, C1QC, IL1B 
       MUC4, LGALS9, CCL2, CXCL16, SOD2, CDH1, CD177, MAFB, LAMC2, FABP2 
PC_ 5 
Positive:  CCL5, CD3E, CXCR4, PTPRC, KLRB1, CD3D, NKG7, HLA-B, CD2, CD8A 
       TRBC1, IFITM1, XCL1, MMP2, GZMK, JAK1, GZMA, GATA3, COL1A1, PDGFRA 
       HLA-C, IFNG, CD3G, CD8B, ICOS, ACTA2, ETS1, GNLY, CD5, ZAP70 
Negative:  MZB1, DERL3, CD79A, CCL3, TNFRSF17, FOS, JUN, CXCL2, HLA-DRA, POU2AF1 
       CCR10, HLA-DPA1, CD27, IRF4, NFKBIA, HLA-DQA1, TNFRSF13B, CXCL1, CDKN1A, HLA-DPB1 
       FCRL5, HLA-DMA, DUSP1, HLA-DRB1, ICAM1, CXCL8, CCL4, ATF3, EGR1, KLF2 

ref <- RunUMAP(ref, dims=1:20)

16:20:37 UMAP embedding parameters a = 0.9922 b = 1.112
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by ‘BiocGenerics’
16:20:37 Read 5725 rows and found 20 numeric columns
16:20:37 Using Annoy for neighbor search, n_neighbors = 30
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by ‘BiocGenerics’
16:20:37 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
16:20:38 Writing NN index file to temp file /project/.tmpRsessions/RtmpQlyPKi/file19fb8dea14db4
16:20:38 Searching Annoy index using 1 thread, search_k = 3000
16:20:41 Annoy recall = 100%
16:20:42 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
16:20:45 Initializing from normalized Laplacian + noise (using RSpectra)
16:20:46 Commencing optimization for 500 epochs, with 235332 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
16:20:55 Optimization finished

DimPlot(ref, label=T, repel=T)